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Abstract 

Model calculations have been performed on the spike-train response of a pair of Hodgkin- 
Huxley (HH) neurons coupled by recurrent excitatory-excitatory couplings with time de- 
lay. The coupled, excitable HH neurons are assumed to receive the two kinds of spike-train 
inputs: the transient input consisting of M impulses for the finite duration (M: integer) 
and the sequential input with the constant interspike interval (ISI). The distribution of 
the output ISI T Q shows a rich of variety depending on the coupling strength and the 
time delay. The comparison is made between the dependence of the output ISI for the 
transient inputs and that for the sequential inputs. 
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1 Introduction 

Neurons communicate by producing sequence of action potentials or spikes. It has been 
widely believed that information is encoded in the average rate of firings, the number 
of action potentials over some suitable intervals. This firing rate hypothesis was first 
proposed by Andrian [|T| from a study of frog, in which the firing rate monotonically 
increases with an increase of the stimulus strength. By applying the firing rate hypothesis, 
the properties of many types of neurons in brain have been investigated and the theoretical 
models have been developed [^]. 

When all action potentials are taken to be identical and only the times of firing of a 
given neuron are considered, we obtain a discrete series of times, {t n }, which is expected 
to contain the information. In the rate coding, only the average of the rate of the interspike 
interval (ISI) is taken into account, and then some or most of this information is neglected. 

In recent years, the alternative temporal coding, in which detailed spike timing is taken 
to play an important role, is supported by experiments in a variety of biological systems: 
sonar processing of bats |J, sound localization of owls [|J, electrosensation in electric 
fish H, visual processing of cats |J@, monkeys [|J and human ||. It is now primarily 
important to understand what kind of code is employed in biological systems: rate code, 
temporal code or others |IJ Q . 



Neural functions are performed in the activity of neurons. Since the Hodgkin-Huxley 



(HH) model was proposed to account for the squid giant axon [O, its property has 



been intensively investigated. Its responses to applied dc [|13|]-|]T7| and sinusoidal currents 



18| [H| have been studied. The HH-type models have been widely employed for a study 



on activities of transducer neurons such as motor and relay neurons, which transform 
amplitude-modulated inputs to spike-train outputs. Regarding the single HH neuron as a 
data-processing neuron, the present author []20| (referred to as I hereafter) has investigated 
its response to the spike-train inputs whose ISIs are modulated by deterministic, semi- 
deterministic (chaotic) and stochastic signals. 

Several investigations have been reported on the property of a pair of the HH neurons 
21|]-|3(J. In the network of two HH oscillators coupled by excitatory couplings without 



time delay, the unit fires periodically in the synchronized state. It is, however, not the 
case when the excitatory couplings have some time delay, for which the anti-phase state 
becomes more stable than the synchronized state ||22|| . Rather, inhibitory couplings with 
substantial time delay lead to the in-phase synchronized states in the coupled HH oscil- 
lators p2| . The similar conclusion is obtained also in the coupled integrate-and-fire (IF) 
oscillators [22,29,31-34]. The phase diagrams for the synchronized state and various clus- 
ter states in the coupled HH oscillators are obtained as functions of the synapse strength 
and the time delay [23-26]. Recurrent loops involving two or more neurons with exci- 
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tatory and/or inhibitory synapses are found in biological systems such as hippocampus 
|35| , neo-cortex |36| and thalamus [|57|| . It is important to make a detailed study on the 



coupled HH neurons, which is the simplest but meaningful network unit. 

In recent years, much attention has been paid to the delayed-feedback systems de- 
scribed by the delay-differential equation (DDE) [38-43]. Their property has been inves- 
tigated with the use of various functional forms for the delay-feedback term in DDE. The 



exposed properties include the odd-harmonic solutions |38flf39|], the bifurcation leading to 
chaos [39-41], the multistability f3"5| , and the chaotic itinerancy ||I2]| ||43| . Among them the 
multistability is intrigue because it may be one of conceivable mechanisms for memory 



storage in biological neural networks. It has been shown by Ikeda and Matsumoto|3T) 
that when the delay time is larger than the response time in the delayed feedback sys- 
tem, information may be stored in temporal patterns. Actually, Foss, Longtin, Mensour 
and Milton [g7| demonstrate this ability in the coupled HH (and IF) neurons with the 
time-delayed feedback. 

The response of DDE is usually discussed by applying the sequential sinusoidal or 
spike-train inputs to non-linear systems. In real neural systems, however, it is not so 
often for neurons to receive such sequential, continuous inputs. Rather, it is expected 
to be more realistic that neurons receive clustered inputs including information to be 
processed. The purpose of the present paper is to investigate the response of the coupled, 
excitable HH neurons to both the transient and sequential spike-train inputs. We adopt 
the recurrent excitatory-excitatory (E-E) couplings between a pair of HH neurons, to 
which we apply the transient spike-train inputs consisting of clustered M impulses for the 
finite duration (M: integer) as well as the sequential inputs with the constant ISI. 

Our paper is organized as follows: In the next §2, we describe a simple neural system 
consisting of neurons, axons, synapses and dendrites, which is adopted for our numerical 
calculation. We present the calculated results in §3: the response of the coupled HH 
neurons to the transient, clustered impulses is discussed in §3.1 and that to the sequential 
spike-train input in §3.2. The dependence of the distribution of the output ISIs on the 
coupling strength and the time delay are studied. The final §4 is devoted to conclusion 
and discussion. 



2 Adopted Model 

We adopt a simple neural system consisting of a pair of neurons which is numbered 1 
and 2. The neurons which are described by the HH model with identical parameters, are 
coupled with the time delay of Tjk (j, k — 1,2) for an impulse propagating from the neuron 
k to the neuron j. This delay time is the sum of conduction times through the axon and 
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dendrite. It has been reported that real biological synapses exhibit temporal dynamics 
of depression or potentiation during neuronal computation [J3|[|3]. We, however, treat 
the synapse as a static unit for a simplification of our calculation. The synapse with 
the coupling strength Cjf. is excitatory, and it is assumed to be described by the alpha 
function [eq. (7)]. 

Dynamics of the membrane potential Vj of the coupled HH neuron j (=1, 2) is described 
by the non-linear DDEs given by 

CdVj(t)/dt = -ipty^m^h^nj) + If* 

+ir t ({v k (t-r jk )}), (i) 



where C = 1 /iF/cm 2 is the capacity of the membrane. The first term of eq. (1) expresses 
the ion current given by 

Ij° n {Vj, rrij, hj, rij) = g Na hj (Vj - V Na ) 

+g K n) (Vj - V K ) + g L (Vj - V L ). (2) 

Here the maximum values of conductivities of Na and K channels and leakage are g-^ a = 
120 mS/cm 2 , g^ = 36 mS/cm 2 and g^ = 0.3 mS/cm 2 , respectively; the respective reversal 
potentials are Vw a = 50 mV, Vk = — 77 mV and Vl = —54.5 mV. The gating variables of 
Na and K channels, rrij, hj and rij, are described by 

drrij/dt = -(a mj + b mj ) rrij + a mj , (3) 

dhj/dt = -(a hj + b hj ) hj + a hj , (4) 

drij/dt = -(a nj + b nj ) rij + a nj . (5) 

The coefficients of a m j and b m j etc. are expressed in terms of Vj (their explicit expressions 
having been given in refs. anc ^ then the variables Vj, rrij, hj and rij are coupled. 

The second term in eq. (1) denotes the external input currents given by 

If = I Bj + A, Strait -tin), (6) 

n 

with the alpha function a(t) given by 

a(t) = (t/r s ) e-*'* Q(t). (7) 

The first term (J SJ ) in eq. (6) is the dc current which determines whether the neuron is 
excitable or periodically oscillating. Its second term expresses the postsynaptic current 
which is induced by the presynaptic spike-train input applied to the neuron 1, given by 

t/iW = KE^-U ( 8 ) 
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In eqs. (2.6)-(2.8), Q(t) = 1 for x > and for x < 0; A s = g s (V a - V s ), g s and V s 
stand for the conductance and reversal potential, respectively, of the synapse; r s is the 
time constant relevant to the synapse conduction, which is assumed to be r s = 2 msec; 
ti n is the n-th firing time of the spike-train inputs defined recurrently by 



Hn+l 



tin T[ n (t[ n ^j , (9) 



where the input ISI T- m is generally a function of t in . For the constant input ISI of T in = 7], 
ti n is given by t- in = nT{ for an integer n. 

When the membrane potential of the j-th neuron Vj(t) oscillates, it yields the spike- 
train output, which may be expressed by 

U oj {t) = K E S{t-t ojm ), (10) 

m 

in a similar form to eq. (8), t j m being the m-th firing time of the neuron j when Vj(t) 
crosses V z = 0.0 mV from below. The output ISI is given by 

T'ojm tojm+1 tojm- (^) 

The third term in eq. (1) which expresses the interaction between the two neurons, is 
assumed to be given by 

if\{v k {t - tjh)}) = E E °jk «(* - T Jk - t okm ). (12) 

Krt) m 

We assume the recurrent excitatory-excitatory couplings with positive Cy given by | 
C21 | = | C12 |= cA s and r 2 i = t 12 = r d . 

As for the functional form of the coupling term of Jj nt ({Vfc(i — Tjk)}), Foss, Longtin, 
Mensour and Milton [27 adopt a simpler form given by 

lf({V k (t - r jk )}) = E Hk V k (t - r jk ), (13) 
K+i) 

taking no account of the synapse, where fij k is the coefficient of the synaptic coupling. 
They discuss the memory storage of the pattern in output spike trains, injecting the input 
information by the initial function, V(t) for t e [—t&, 0), whereas in our calculation input 
information is given by 7| xt [eq. (6)]. 

Differential equations given by eqs. (l)-(5) including the external current and couplings 
given by eqs. (6)-(12) are solved by the forth-order Runge-Kutta method. The calculation 
for each set of parameters is performed for 2 sec (200,000 steps) with the integration time 
step of 0.01 msec with double precision. The initial conditions for the variables are given 
by 

Vj(t) = -65 mV, mj (t) = 0.0526, hj(t) = 0.600, 

rij{t) = 0.313, for j = 1, 2 at t = 0, (14) 
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which are the rest-state solution of a single HH neuron (Cjk = 0). The initial function for 
Vj(t), whose setting is indispensable for the delay-differential equation, is given by 

Vj(t) = -65 mV for j = 1, 2 at t e [-r d , 0). (15) 

For an analysis of asymptotic solutions, we discard results of initial 1000 msec (100,000 
steps). 



3 Calculated Resuls 

In the present study, we consider only the excitable HH neurons by setting J s , = and 



2 r n r o / 2 



A = g s (V a - V 3 ) = 40/iA/cni for g s = 0.5mS/cni , V a = 30 and V s = -50 mV |g. Our 
model has additional three parameters, 7], t<± and c. We treat them as free parameters to 
be changed because the values of ISI and the time delay observed in biological systems 



distribute in a fairly wide range [46| 



3.1 Transient Spike- Train Inputs 

Let us first investigate the response to the transient, clustered spike-train inputs consisting 
of M impulses. We have studied in I, the transient response of a single HH neuron to 
spike-train inputs consisting of M = 2 — 5 impulses with T; = 5, 10 and 20 msec (see Fig. 
20 of I). In the case of 7] = 20 msec, we get T Q = 20 msec and the number of output 
pulses is the same as that of input pulses. On the contrary, in the cases of TJ = 5 and 10 
msec, the ISI of output is generally larger than that of input because of its character of 
the low-pass filter, and the number of output pulse is not necessary the same as that of 
input pulse. 

Figure 1 shows the example of the time courses of input (U{), output pulses (U j), the 
total postsynaptic current (Ij = Ij xt + Ij nt ) and the membrane potential (Vj) with M = 3, 
T\ = 20, Ta = 10 msec and c = 1.0 for the E-E coupling (c > 0). The first external pulse 
applied at t = yields the firing of the neuron 1 after the intrinsic delay of r ;i ~ 2 msec. 
The emitted impulse propagates the axon and reaches the synapse of the neuron 2 after 
T21 = 10 msec. After a more delay of an intrinsic 712 ~ 2 msec, the neuron 2 makes the 
firing which yields the input current to the neuron 1 after a delay of r 12 = 10 msec. The 
input pulses trigger the continuous oscillation in the coupled HH neurons with the output 
ISI of T Q = 24.10 msec. The time dependence of the output ISI of the neuron 1 and 2 
are plotted by solid and dashed lines in Fig. 2(a), respectively. We note that T ol and T o2 
start from the values of 20.00 and 19.96 msec, respectively, and soon become the value of 
24.10 msec. 
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Figures 2(b) and 2(c) show similar plots for different values of r d = 13.75 and 20 msec, 
respectively. In the case of r d =13.75 msec, output ISIs start the oscillation with T o =20.00 
and 12.36 msec and asymptotically approach the value of 15.93 msec. On the contrary, 
in the case of Td=20 msec, the oscillation of T G starting at t — continues with the 
asymptotic values of T = 18.51 and 25.59 msec. In the following subsections, we will 
discuss the dependence of output ISIs on the time delay and the coupling strength. Since 
its behavior of the spike-train outputs of the neurons 1 and 2 is similar, we hereafter take 
into account only that of the neuron 1 otherwise noticed. 



3.1.1 The time-delay dependence 

Now we study how the output ISIs are determined. When the coupling strength is suf- 
ficiently strong for inputs to trigger output impulses and when the feedback time 7ft, is 
larger than the duration of clustered impulses (i.e. Tr, = 2r d + m + r i2 > (M — 1) Tj), we 
get two values of T G given by 

7^(1) _ 7"! 

T Q (2) = Tfb — (M — 1) 7] 

= 2r d + r il + r i2 -(M-l)T i . (16) 

On the other hand, when the feedback time is shorter than input-pulse duration (2r d + 
Tii + Tis < (M-l)71), we get 

T« = 71 0(M-3), 
T Q (2) =| £Tfb - mT; | 

=| t (2r d + Tii + r i2 ) -mTl |, (17) 

where integers i and m satisfy 1 < £ < [(M - 1)71/7^] + 1 and < m < M - 1, [ • ] is 
the Gauss sign and is vanishing for M < 2. 

Figures 3(a) and 3(b) show the calculated time-delay dependence of T for c = 1.0 
and 1.6, respectively. Filled and open circles denote T Q s of the transient (t < 1000 msec) 
and asymptotic solutions (t > 1000 msec), respectively. As was shown in Figs. 2(a)-2(c), 
the output ISIs of the asymptotic solutions are T o =24.10 for r d = 10 msec, T Q =15.93 for 
r d =13.75 msec, and T Q =18.51 and 25.59 msec for r d =20 msec. We note that the behavior 
of T Q strongly depends on the value of r d . In order to see their detailed structures, we show 
in Figs. 4(a)-4(c), enlarged plots of the narrow regions in Fig. 3(b). Figures 3(a) and 3(b) 
show three main branches expressed by T Q ~ 2r d + 5, T Q ~ 2r d — 15 and T Q ~ 2r d — 35, 
which are obtainable for a pair of integers of (£, m) =(1,0), (1,1) and (1,2), respectively, 
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by eq. (17) with t x \ = = 2.5 and 7] = 20 msec. Figure 4(a), in which the narrow region 
of 10 < r d < 20 msec in Fig. 3(b) is enlarged, shows an additional branch of a single ISI 
given by T ~ r d + 2 at 11.9 < r d < 12.3 and 13.7 < r d ~ 19 msec. Furthermore we note 
in Fig. 4(b) which shows an enlarged plot at 20 < r d < 30 msec in Fig. 3(b), a branch of 
multiple ISIs given by T Q ~ 0.5r d + 5 at 21.9 < r d < 28.5 msec. These r d dependences of 
T Q cannot be explained by eq. (17), and may be harmonics of the fundamental ISI with 
the period of 2r d . The r d dependence of T Q for c = 1.6 shown in Fig. 3(b) is similar to 
that for c = 1.0 shown in Fig. 3(a), except an additional branch given by T Q ~ 0.5r d + 4 at 
12.3 < r d < 13.7 msec. The narrow region of 16 < r d < 20 msec in Fig. 4(a) is enlarged 
in Fig. 4(c), where the ISI of the asymptotic solution shows the stair-like structure. 



3.1.2 The coupling-strength dependence 

Figures 5(a) and 5(b) show the c dependence of the output ISI of the transient (filled 
circles) and asymptotic solutions (open circles). As was shown in Figs. 2(a) and 2(b), the 
ISI of the asymptotic solutions with c = 1.0 is T Q =24 msec for r d = 10 msec and T Q =15.93 
msec for r d = 13.75 msec. Figure 5(a) shows that as the coupling strength becomes weak, 
T Q is increased because of the integrator character of the HH neuron. A similar effect is 
obtained also in a single HH neuron, in which the output ISI becomes larger for smaller 



spike-train inputs [20|. Figure 5(b) shows that ISIs for the transient solutions fluctuate 
around that for the asymptotic solution as expected. The enlarged plot for 1.5 < c < 1.7 
of Fig. 5(b) is given in Fig. 6, where a discontinuous change in T Q is clearly realized at 
c = 1.61 msec. For c < 0.2 neurons emit only three impulses, returning to rest without 
oscillations. 



3.2 Sequential Spike- Train Inputs 

Next we discuss the response to the sequential spike-train. Our calculations in I show that 
when an isolated HH neuron (c = 0) receives the sequential inputs with the constant ISI 
of 7], it behaves as a low-pass filter: it emits the spike train with T Q > 10 msec for 7] < 12 
msec while for 7] > 12 msec its output ISI is given by T Q = 7j (see Fig. 7 of ref. 20). 
This response may be modified when the coupling is introduced to a pair of HH neurons. 
Figure 7 shows the time courses of input (Ui), output (U Q j), the total postsynaptic current 
(Ij = 7J xt + 7j nt ) and the membrane potential (Vj) for 7] = 20 msec, r d = 10 msec and 
c = 1.0, which are the same parameters adopted for the clustered inputs shown in Fig. 
1. The output ISI in Fig. 1 is 24.1 msec while that in Fig. 7 is 20 msec which is the 
entrained value with input ISI. The response behavior of the coupled neurons strongly 
depends on the parameters of c, r d and T\. 
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3.2.1 The time-delay dependence 

Figures 8(a) and 8(b) show the r d dependence of the distribution of T Q for c = 1.0 and 



1.6, respectively, in the asymptotic solution of the sequential inputs [|7|. The calculations 
in Figs. 8(a) and 8(b) are performed with the same parameters of Tj and c in Figs. 3(a) 
and 3(b), respectively. The r d dependence of T Q in Fig. 8(a) [Fig. 8(b)] is quite different 
from that in Fig. 3(a) [Fig. 3(b)]. Figures 8(a) and 8(b) have the bifurcation structure, 
as commonly observed in systems with the delayed feedback [Q. In order to see more 
the detailed structure of the bifurcation, we show, in Fig. 9, the enlarged plot for the 
range of 21 < t& < 26 msec between the dotted, vertical lines in Fig. 8(a). The region 
sandwiched by verical dotted lines in Fig. 9(a) (21 < r<j < 26 msec) is further enlarged 
in Fig. 9(b). Figures 9(a) and 9(b) clearly show the bifurcation as changing The Td 
dependence for c = 1.6 shown in Fig. 8(b) is similar to that for c = 1.0 shown in Fig. 
8(a), and its enlarged plot also exhibits the bifurcation (not shown). 



3.2.2 The coupling-strength dependence 

The calculated c dependence of the distribution of T Q with = 10 and 13.75 msec for 



T\ = 20 msec are shown in Figs. 10(a) and 10(b), respectively [47]. The adopted values 
of Ti and t& in Figs. 10(a) and 10(b) are the same as those in Figs. 5(a) and 5(b), 
respectively. The output ISI for the sequential inputs shown in Fig. 10(a) is 20 msec 
(= Ti) independent of the coupling constant, while that for the clustered inputs shown in 
Fig. 5(a) decreases monotonically as the c value is decreased. We note in Fig. 10(b) that, 
as increasing the c value, the distribution of the output ISIs for r = 13.75 msec exhibits 
the bifurcation. In order to investigate the phenomenon in more detail, we show, in Fig. 
11, the enlarged plot for the range of 0.6 < c < 1.2 sandwiched by the dotted, vertical 
lines in Fig. 10(b). 

A cycle whose output ISIs almost continuously distribute, is expected to be chaotic 
although in the strict sense, the distribution of our T Q s never becomes continuous because 
they are quantizedby the integration time step of 0.01 msec. Among many candidates of 
chaos-like behavior in Figs. 11, we pay our attention to the result of c = 0.95, for which 
the Lorentz plot (return map) of its T Q is shown in Fig. 12(a) (calculations are performed 
for 20 sec of two million steps). The output ISIs seem to distribute on the folded ring. 
When these points are connected by lines in the temporal order, the inside of the ring is 
nearly filled by them. In order to examine the property of this cycle, we calculate the 
correlation dimension v given by [^8 

v = hm 



+ o log e 



10 



with 

N 

C(e) = N- 2 ©(e- \X m -X n |), (19) 

m,n=l 

X m (T om , T om +i, , 7om+fc— l)j (20) 

where C(e) is the correlation integral, X m is the /c-dimensional vector generated by T om , 
N the size of data, and ©(•) the Heaviside function. Figure 12(b) shows the logC(e)-loge 
plot for various embedding dimensions k calculated for the cycle shown in Fig. 12(a) with 
the data size of iV ~ 1200. We note that C(e) behaves as C(e) oc e v with the correlation 
dimension of v — 0.94 ±0.02 for small e (0.01 = e -4 6 < e < e°). When the relevant spike- 
train output given by U i(t) [eq. (10)] is Fourier transformed, it spectrum shows a broad 
distribution. These suggest that the cycle shown in Fig. 12(a) may be chaos, although we 
cannot draw any definite conclusion until a detailed calculation of its Lyapunov spectrum 
is performed, related discussion being given in §4. 



4 Conclusion and Discussion 

We have performed model calculations of the spike-train responses of a pair of coupled 
HH neurons, applying the two types of inputs of the transient and sequential spike-train 
impulses. Calculations for the transient inputs shown in Figs. 3(a) and 3(b) [Fig. 5(a) and 
5(b)] are performed with the same model parameters as those for the sequential inputs 
shown in Figs. 8(a) and 8(b) [10(a) and 10(b)]. When we make a comparison of the 
response to transient inputs with the corresponding result to sequential inputs, we notice 
the difference and similarity between them. When we regard a neuron as a data-processing 
element, the relation between input and output ISIs is one of the important factors for 
its quality. Our previous calculation in I shows that a single HH neuron emits a single 
ISI of T Q = 7] for 7] < 12 msec whereas for shorter ISI of 7] < 12 msec it emits multiple 
ISIs of T Q > 10 msec (see Fig. 7 of I). Figures 13(a) and 13(b) show the 7] dependence 
of T Q of coupled HH neurons for the transient and sequential inputs, respectively, with 
r d = 50 msec and c = 1.0. Dashed lines in Figs. 13(a) and 13(b) are obtained with the 
use of eq. (17) for a pair of integers (£,m) shown in the brackets. It is apparent that the 
distribution of T Q in Fig. 13(a) is not the same as that in Fig. 13(b), but they are partly 
similar. 

On the theoretical point of view, the sequential input is taken as the limit of M — > oo 
of the M-impulse clustered input. In order to understand the transition of the response 
behavior as increasing M, we plot, in Fig. 14, the time dependence of T Q for this set of 
parameters by changing the M value. For M = 3, T Q oscillates with the values of 19.6, 
20.0 and 64.54 msec, as mentioned before. The calculated T Q for M=4 are 18.7, 19.8, 
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20.0, and 45.6 msec, and those for M = 5 are 19.9, 20.0 and 24.2 msec. For M = 10, T D 
remains 20 msec until t ~ 200 msec, after which T Q oscillates with the values of 19.9, 20.0 
and 24.2 msec. In the limit of M — > oo corresponding to the sequential inputs, the state 
with T Q = 20 msec continues from t = to oo. Thus as increasing M, the time region 
of T Q = 20 msec is increased. Figure 15 shows the similar plot of the time dependence 
of the distribution of T D for various M with Tj = 20, = 13.75 and c = 0.95, for which 
the sequential input leads to the chaotic behavior, as was discussed in Sec. 3.2 (see Fig. 
12(a)). In the case of M = 3, we get the oscillation in T Q which asymptotically approaches 
the value of 15.97 msec. In the case of M — 10 (50), the chaotic behavior is realized at 

< t ~ 180 (0 < t ~ 980) msec during the application of inputs. After inputs are 
switched off, the output ISI gradually approaches the asymptotic value of 15.97 msec. In 
the limit of M — > oo, the chaotic oscillation eternally continues. 

We have shown in §3.2 that the cycle of the output ISIs shown in Fig. 12(a) may be 
chaos because of its correlation dimension of v ~ 0.94 derived from the logC(e)-loge plot 
in Fig. 12(b). This is not surprising because the response of single HH neurons to some 
kinds of external inputs may be chaotic [18l ||191||20|| . In particular, it has been shown in 

1 that the response of a single HH neuron may be chaos when the ISI of the spike-train 



input is modulated by the sinusoidal signal: f20[ 

Ti(t) = d + di sin(27rt/T p ), (21) 

where do denotes the average of 71 (t), d\ the magnitude of the sinusoidal modulation, and 
T p its period. Figure 16(a) shows the Lorentz plot of the output ISIs of the single HH 
neuron (c = 0.0) receiving sequential inputs modulated by sinusoidal ISIs [eq. (21)] with 
d = 2di = 20 and T p = 100 msec (see Fig. 9(d) of I, where points in the Lorentz plot are 
connected by lines in the temporal order). We note that T Q s distribute on the deformed 
ring. From the log C (e)-log e plot (not shown) of this cycle, we get its correlation dimension 
of v <~ 1.04. We apply this sinusoidal spike-train input to the coupled HH neurons with 
Td = 10 msec and c = 1.0, whose Lorentz plot is shown in Fig. 16(b). Its structure is rather 
different from that shown in Fig. 16(a). Actually the correlation dimension of this cycle 
for the coupled HH neurons is v ~ 1.83, which is different from and larger than v ~ 1.04 
of the cycle shown in Fig. 16(a) for the single HH neuron. From similar calculations for 
the coupled HH neurons, we obtain the correlation dimensions of v ~ 0.95 for = 5 msec 
and c = 1.0, and v ~ 1.03 for Td = 10 msec and c = 0.5. These results clearly show that 
the correlation dimension of the output ISIs depend not only on the model parameters (c 
and Td) of the coupled HH neurons but also on v u the correlation dimension of input ISIs 
{i>i = for the constant ISI and V{ — 1 for the sinusoidally modulated ISI). We expect 
that spike-train inputs with larger v\ lead to spike-train outputs with larger v. One of 
the disadvantages of the present calculation of the correlation dimension is a lack of the 
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data size of iV ~ 1200 with million-step calculations. A more accurate analysis requires 
a larger size of data and then a computer with the larger memory storage. 

Next we discuss the time correlation r 12 (r) between the membrane potentials, V\ and 
V 2 , of the neurons 1 and 2, defined by 



ri 2 (r) = / [Vxit)- < V 1 (t) >] [V 2 (t + r)- < V 2 (t) >] dt, (22) 



where the bracket denotes the time average, and t a = 1000 and tb = 2000 msec are 
adopted for our calculation. Figure 17(a) shows the result for the case of the sequential 
input to the coupled HH neurons with T ; = 20, = 10 msec and c = 1.0 (see Fig. 7 for 
the time courses of Vy and V 2 ). In this case we obtain the constant T Q = 20 msec as was 
discussed in § 3.2, and then r 12 (r) shown in Fig. 17(a) has peaks at r = 12.04 + 20 n 
msec (n: integer) with the period of 20 msec, as expected. We are interested in the time 
correlation for the case when the distribution of T Q is chaotic. Results for such cases are 
shown in Figs. 17(b) and 17(c). We have discussed in § 3.2 that the cycle of T Q depicted 
in Fig. 12 may be chaotic. Figure 17(b) shows the result of this case for the E-E coupling 
with T; = 20, Td = 13.75 msec and c = 0.95. We note that Tu(t) has peaks at r=-45.34, 
-30.69, -16.13, 0.0, 16.07, 30.54, 48.17,... msec with the period of about 16 msec, which 
is the sum of r d and m. More evident peaks are found in Fig. 17(c) showing also the 
chaotic case discussed in the preceding paragraph: the E-E coupled neurons receiving the 
sinusoidal inputs given by eq. (21) with d = 2d\ = 20, T p = 100, T| = 20, r d = 10 
msec and c = 1.0 [see Fig. 16(b)]. We note peaks in Ti 2 (r) at r=-37.63, -25.03, -12.52 
0.0, 12.72, 25.30, 37.88, ... msec with the period of about 12.6 msec. These results are 
not modified even when the initial condition of one of the HH neurons is slightly changed 
from the values in eq. (14). It is interesting that the synchronization is well preserved 
between the coupled HH neurons in the chaotic state |^9| . 



A fairly large variability (c v = 0.5 ~ 1.0) has been reported for spike trains of non- 
bursting cortical neurons in VI and MT of monkey |50| . It is possible that when the 
appreciable variability in neuronal signals is taken into account in our calculations, much 
of the fine structures in the c— and rj-dependent distributions of T Q will be washed out. In 
order to study this speculation, we apply the spike-train input with ISI whose distribution 
is given by the gamma distribution defined by 



P(T) = s r T e~ si / r(r), (23) 

where T (r) is the gamma function. The average of input ISI is given by fi\ = r/s, its 
root-mean-square (RMS) by o\ = y/r/s and its variability by c vl = o\j [i\ = 1/v 7 ^- Figure 
18 shows the ra dependence of the mean (/i ) and RMS values (er ) of the output ISIs for 
c v ; = 0.0 (dashed curves) and c v i = 0.43 (solid curves) with T\ = 20 msec and c = 1.0. 
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Note that a Q provides us with the measure of the width of the distribution of T Q . The 
distribution for c vl = has a fine structure reflecting the strong Td dependence of T Q [see 
Fig. 8(a)]. This fine structure is, however, washed out for c vi = 0.43, as expected. 

Finally we discuss the relevance of the calculated properties to biological experiments. 
Many experimental data have shown the complex behavior of electoencephalographic 
(EEG) waves in brain. The macroscopic characteristics of their activity are aperiodic 
and unpredicable oscillations with amplitude histograms that are near Gaussian, auto- 
correlation functions that rapidly approach zero and intermittent burst of oscillations 
having spectral peaks It has been reported that the activity of EEG in the olfactory 
bulb shows the significance of chaos in an animal's motivated behaviors ]5^| . The complex 
behavior of EEG is nothing but the reflection of that of action potentials generated by 
neurons. It has been shown that neurons in different regions of the brain have different 
firings property. In hippocampus, for example, gamma oscillation (20 ~ 70 Hz) occurs 
in vivo, following sharp waves |35|]. In neo-cortex, gamma oscillation is observed under 
conditions of sensory signal as well as during sleep |36| . In thalamus burst firings are found 
during slow-sleep, and single spiking is found during arousal |37|]- One of the reasons of 
this variety of firings is that different class of neurons has different ion conductances. 
Physiological experiments have shown that these biological systems include recurrent 
loops connecting two or more neurons with excitatory and/or inhibitory synapses. It is 
conceivable that the distributed processing of brain function may be due to differences not 
only in ion conductances of the neuron but also in synaptic strength and in delay times 
of axons and dendrites connecting neurons. Although many theoretical studies have been 
made, the origin of the complexity in neuron firings has not been well clarified at the 
moment. We should note that synaptic strengths may be modified by Hebb's learning 
rule, which changes the state of the network including given synapses. Our calculations 
for a pair of HH neurons, which is a simplest, plausible model simulating recurrently 
connected network, have demonstrated that depending on the coupling strength and the 
time delay, the coupled HH neurons show a much variety, yielding not only regular spike 
trains but also irregular (chaotic) impulses. We hope that our calculations might have 
some relevance to the complex activities in real, biological systems. 
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Figure Captions 

Fig.l The time dependence of the clustered input (£/;), output (U Q j), the total postsy- 
naptic current (/,-) and the membrane potential (Vf) with M = 3, T; = 20, r d = 10 msec 
and c = 1.0. The result of V 2 is shifted downward by 200 mV and scales for U h U j and 
Ij are arbitrary. 

Fig. 2 The time dependence of the output ISI (T OJ ) of neuron 1 (filled circles) and 2 (open 
circles) for (a) r d = 10, (b) 13.75 and (c) 20 msec for clustered inputs with M = 3, T| = 20 
msec and c = 1.0. 

Fig. 3 The r d dependence of the distribution of T G of (a) c = 1.0 and (b) 1.6 for the 
clustered inputs with M = 3 and 7] = 20 msec. Filled and open circles denote the results 
of the transient (t < 1000 msec) and asymptotic solutions (t > 1000 msec), respectively. 
The dashed lines are expressed by the equations written beside the lines. The enlarged 
plots of the regions between dotted, vertical lines in Fig. 3(b) are shown in Figs. 4(a)-4(c) 
(see text). 

Fig. 4 The enlarged plot of the r d dependence for (a) 10 ~ r d ~ 20 msec, (b) 20 ~ r d ~ 30 
msec, and (c) 16 ~ r d ~ 20 msec for M = 3, 7- = 20 and c = 1.6@ (see Fig. 3(b)). 

Fig. 5 The c dependence of the distribution of T G for (a) r d = 10 and (b) 13.75 msec to the 
clustered inputs with M = 3 and T\ = 20 msec. Filled and open circles denote the results 
of the transient (t < 1000 msec) and asymptotic solutions (t > 1000 msec), respectively. 
The enlarged plot of the regions between dotted, vertical lines in Fig. 5(b) is shown in 
Fig. 6 

Fig. 6 The enlarged plot of the c dependence of the T Q for M = 3, T ; = 20 and r d = 13.75 
msec (see Fig. 5(b)). 

Fig. 7 The time course of sequential input (£/;), output (U Q j), the total postsynaptic 
current (Ij), and the membrane potential (Vj) with 7] = 20, r d = 10 msec and c = 1.0. 
The result of V 2 is shifted downward by 200 mV and scales for U h U j and Ij are arbitrary. 

Fig. 8 The r d dependence of the distribution of T G of (a) c = 1.0 and (b) 1.6 for the 
sequential input with T\ = 20 msec. The enlarged plot of the region between dotted, 
vertical lines in Fig. 8(a) is shown in Fig. 9. 
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Fig. 9 The enlarged plot of the Td dependence of T Q at (a) 21 < Td < 26 msec and (b) 
23 < Td < 24 msec [the region sandwiched by vertical dotted lines in 9(a)] for 7j = 20 
msec and c = 1.0 (see Fig. 8(a)). 

Fig. 10 The c dependence of the distribution of T Q for (a) r d =10.0 and (b) 13.75 msec to 
the sequential inputs with 71=20 msec. The enlarged plot of the regions between dotted, 
vertical lines in 10(b) is shown in Fig. 11. 

Fig. 11 The enlarged plots of the c dependence of T Q for 7j = 20 and r d = 13.75 msec (see 
Fig. 10(b)). The arrow denotes the c value for which the Lorentz plot is shown in Figs. 
12(a). 

Fig. 12 (a) The Lorenz plot of T om for c = 0.95 with 7] = 20 and Td=13.75, the computa- 
tion being performed for 20 sec (two million steps), (b) The correlation integral C(e) of 
the cylce shown in (a) as a function of e in the log-log plot for various dimensions k, the 
dashed line denoting C(e) oc t v with the correlation dimension of v = 0.94 [eqs.(3.3)-(3.5)]. 

Fig. 13 The 7] dependence of the distribution of T G for (a) the clustered input (M = 3) 
and (b) sequential spike-train input with r d = 50 msec and c = 1.0. Filled and open 
circles in (a) denote the results of the transient (t < 1000 msec) and asymptotic solutions 
(t > 1000 msec), respectively, while in (b) filled circles express the results of asymptotic 
solutions (t > 1000 msec). Dashed lines are expressed by a pair of integers of (£, m) in 
eq. (17) (see text). 

Fig. 14 The time dependence of T D for the clustered impulse inputs with M = 3, 10, 50 
and oo with 7] = 20, Td = 50 msec and c = 1.0, results of M=3, 10 and 50 being shifted 
upward by 60, 40 and 20 msec, respectively. The arrows denote the time below which the 
inputs are continuously applied. 

Fig. 15 The time dependence of T Q for the clustered impulse inputs with M = 3, 10, 50 
and oo with 7] = 20, r d = 13.75 msec and c = 0.95, results of M=3, 10 and 50 being 
shifted upward by 30, 20 and 10 msec, respectively. 

Fig. 16 The Lorenz plots of T om of (a) the single HH neuron (c = 0.0) and (b) the coupled 
HH neurons (rd = 10 msec and c = 1.0) receiving spike-train inputs whose ISIs are 
modulated by sinusoidal signal given by eq. (21) with d = 2di = 20 and T p = 100 msec 
(see text). 
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Fig. 17 The time correlation r^fV) between the membrane potentials of the neurons 1 
and 2 [eq. (22)] for (a) the constant-ISI input with 7] = 20, Td = 10 msec and c = 1.0, 
(b) that with 7] = 20. = 13.75 msec and c = 0.95, and (c) the sinusoidal spike-train 
input given by eq. (21) with d\ = 2d 2 = 20, T p = 100, r d = 10 msec and c = 1.0. The 
results of (b) and (c) are shifted downward by 1.0 and 2.0, respectively (see text). 

Fig. 18 The r d dependence of the mean (/i ) and rms (er ) of output ISIs of the coupled 
HH neurons (c = 1) receiving sequential inputs of =< 7] >= 20 msec with c vi = 0.0 
(dashed curves) and c vi = 0.43 (solid curves) (see text). 
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